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Abstract. Consider the overdamped limit for a system of interacting particles in the presence of 
hydrodynamic interactions. For two-body hydrodynamic interactions and one- and two-body poten- 
tials, a Smoluchowski-type evolution equation is rigorously derived for the one-particle distribution 
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1. Introduction. 

1.1. Review of existing dynamic density functional theory. Several prob- 
lems in condensed matter physics such as colloidal suspensions and polymers can often 
be described as systems of interacting Brownian particles, either in the presence or 
absence of hydrodynamic interactions HZ] . Hydrodynamic interactions are due 
to forces on the colloid particles caused by flows in the suspending fluid (referred to 
as the bath), which are generated by the motion of the colloidal particles, and can 
be thought of as generalised friction forces. Such systems of interacting Brownian 
particles can be described either in phase space, when both positions and momenta 
of the particles are taken into account, or in configuration space, when inertial effects 
are neglected and only the position of the Brownian particles is taken into account. 
The evolution of the phase space distribution function is described by the Kramers 
equation [331 134j . On the other hand, the evolution of the distribution function in 
configuration space is governed by the Smoluchowski equation [7TJ [51] . 

The Smoluchowski equation can be derived from the Kramers equation in the 
overdamped, i.e. high friction, limit. In this limit the velocity of the particles ther- 
malizes quickly, i.e. the velocity distribution converges quickly to a Maxwellian and 
the momentum variables can be eliminated through an appropriate adiabatic elimi- 
nation procedure. This procedure is now well understood, both for a single Brownian 
particle as well as for systems of interacting particles, and both in finite as well as 
infinite dimensions pi [13l |45] [57| [30] . 

Whilst the derivation of the Smoluchowski equation was already discussed by 
both Klein and Kramers, their approach was largely heuristic and a rigorous theory 
was not introduced until later [45] . Systematic adiabatic elimination techniques were 
introduced in the '70s and applied to the problem of the rigorous derivation of the 
Smoluchowski equation for the Kramers equation, e.g. by Wilemski [73j and Titu- 
laer [BEIISH]- In particular, Titulaer considered the fully-interacting iV-body linear 
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Kramers equation and used a multiplc-time-scale or Chapman-Enskog expansion in 
the friction constant to systematically derive the iV-body Smoluchowski equation and 
its corrections. These systematic adiabatic elimination procedures can be understood 
in the context of singular perturbation theory for Markov processes |48l or, more gen- 
erally, in the framework of multiscale methods (e.g. 50 ). A pedagogical discussion 
of the multiple-time-scale technique was given by Bocquet [8], for the case of non- 
interacting particles and for the one-body reduced distribution function. See also [501 
Ch. 11] and [121 Ch. 8]. It is worth noting that in all these derivations, the equations 
for the one-body distribution function are linear. 

Both the Kramers and Smoluchowski equations describe the full iV-body dynam- 
ics, and, although they are linear, they are not well-suited to computation. This is 
due to the large number of variables, and whilst the derivation of the Smoluchowski 
equation from the Kramers equation is of fundamental interest, the reduction from 
67V to 3JV variables by eliminating the momentum variables is insignificant in terms 
of computational complexity. However, a further simplification arises by integrating 
out over the positions (and momenta) of all but one particle, which then allows to 



obtain the dynamics of the reduced distribution functions 0153]; see also (2.2 1. One 
of the main goals of statistical mechanics and kinetic theory is the derivation of closed 
equations for these reduced distribution functions, and in particular for the one-body 
distribution function. If there are no inter-particle interactions, i.e. the only force 
comes from an external potential, and hydrodynamic interactions between particles 
are neglected, this reduction procedure results in one-particle versions of the Kramers 
and Smoluchowski equations. In the more general case where such inter-particle effects 
may not be neglected, the equations must be closed by choosing a suitable approxi- 
mation of the higher-body densities in terms of the one-body distribution function, 
e.g. a mean field approximation. Such a description is the ultimate aim of dynamic 
density functional theory (DDFT). 

Consider a system of N interacting particles with iV-body distribution func- 
tion /W(ri,pi, . . . ,rjv,pjv,r), which gives the probability of finding particles at 
ri, . . . ,rjy with momenta pi, . . . , pjv at time r. The derivation of a self-consistent 
DDFT requires expressing the full ./V-body distri butio n function f( N ^ in terms of the 



one-body reduced distribution f^\ri, pi, r) (sec < 2.2 )) or, if starting from the Smolu- 
chowski equation, in terms of p(ri, r) := j dpi/^^ri, pi, r). Whilst it is known that 
f( N ) (and thus all properties of the system including all lower rt-body distributions) 
is given by a unique functional of p, both in 32, 43.! and out |14j of equilibrium, in the 
general case this functional is unknown. However, much work has been done for the 
equilibrium case (density functional theory, or DFT which allows for an accurate de- 
scription of the microscopic properties of a fluid in terms of its density distribution; see 
[331 [SH] for early work and e.g. (7H[7S] for recent overviews), and there exist accurate 
functionals e.g. Rosenfeld's fundamental measure theory for hard spheres [59l [60] [61] 
and the mean field approximation \37\ [1], mentioned earlier, which becomes exact 
for soft interactions at high densities. DFT represents one of the most widely used 
methods in condensed matter physics for the study of the microscopic structure of 
non-homogeneous fluids within the framework of equilibrium statistical mechanics. It 
offers an increasingly popular compromise between computationally costly molecular 
dynamics simulations and various phenomenological approaches. It has been used 
to describe a wide variety of physical settings, ranging from polymers |36j . liquid 
crystals [T5J and molecular self-assembly [57J 12H] to interfacial phenomena including 
wetting transitions on substrates [551 ES] • 
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We now discuss how one moves from DFT to DDFT, where the system lies away 
from equilibrium. DDFT is also a popular approach in condensed matter physics, 
and has been applied to a wide range of problems including spherical colloids without 
hydrodynamic interactions in both configuration [3U GM HOI HI] and phase space 
[38l [3] , dense atomic liquids [2] , anisotropic colloids [56] , and inhomogeneous granular 
fluids [41]. The effects of inertia [42] and hydrodynamic interactions [Ml [52] have 
also been studied. However, none of the formalisms derived so far is rigorous. The 
relationships between different approaches are summarized in Figure [j~T| We note, 
in particular, that the two routes to obtain the one-body Smoluchowski equation 
give, in general, different formulations. These derivations can be divided into four 
cases, firstly by whether they start from the Kramers or Smoluchowski equation, and 
secondly by whether or not they include hydrodynamic interactions. 

1.2. Starting from the Kramers or the Smoluchowski equation. As men- 
tioned earlier, it is expected that the Smoluchowski equation is valid in the over- 
damped limit, whereas for intermediate and small values of the friction coefficient the 
Kramers equation should be used. Perhaps the most common additional approxima- 
tion is to ignore the effects of the hydrodynamic interactions between the particles 
(for exceptions, see [531 Whilst this may be acceptable in a very dilute system, 
such interactions decay only polynomially slowly with inter-particle distance, and are 
thus long-range and important in many applications |20j . 

When starting from the Smoluchowski equation and neglecting hydrodynamic 
interactions, it suffices to employ the adiabatic approximation, first introduced by 
Marconi and Tarazona [39l [40] . At equilibrium, Mcrmin's proof [43] shows that there 
exists a unique functional of p, J-" ox [/o], called the excess free energy functional, which 
exactly determines the contributions from the many-body potentials (which a priori 
involve higher-order reduced distributions). It then remains to determine accurate, 
generally empirical, approximations to the unknown functional •7 r G x[/>]- The adiabatic 
approximation assumes that the same relationship holds away from equilibrium. This 
is equivalent to assuming that the non-equilibrium n-body distributions are identi- 
cal to those in an equilibrium system with the same instantaneous density p. This 
approximation has proven accurate in a range of systems [TJ [H 1551 162j . 

If hydrodynamic interactions are included, this approximation is insufficient. This 
is because there are no hydrodynamic effects at equilibrium. Instead, at least for 2- 
body interactions, one uses the identity (r, r', r) = p(r, t)p(r', r)g(r, r'; [p]), where 
g is a pair-distribution function, while the function g — 1 is known as the pair correla- 
tion function (it provides a measure of the distance over which particles are correlated; 
for an ideal gas g = 1) and assumes that a good approximation to g is known |54| 
(often g can be approximated with different methodologies, such as the BBGKY hi- 
erarchy or the Ornstein-Zernike equation). Note in particular that g is a functional 
of p. 

When starting from the Kramers equation an additional problem is encountered. 
In this case one obtains an infinite hierarchy of equations for the evolution of the 
momentum moments of /^(r, p, r), i.e. for J dpp^p^p^ 3 (r, p, r) with a,j > 0, 
o,j = n. To obtain closure, one must truncate this hierarchy at a given level, 
which requires the approximation of higher moments. For example, the standard 
truncation at the velocity [n — 1) level (i.e. the same level of description as the 
Navier-Stokes equations), one must control terms (in appropriate units) of the form 
/dp (p ® p — l)/*- 1 -* (r, p, t), where 1 is the 3x3 identity matrix. 

At equilibrium, this term vanishes. However, it is analogous to the kinetic energy 
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Fig. 1.1. Flow diagram of the various approaches used to obtain one-body evolution equations 
and DDFTs from the full underlying dynamics. Arrows indicate the interconnectedness of the differ- 
ent approaches. Thick boxes/arrows: this work. Thin boxes/arrows: previous approaches. Dashed 
boxes/arrows: note that the two routes produce different one-body Smoluchowski equations when 
hydrodynamic interactions are included, although both approaches are accui •ate to 0(e 2 ). Text on 
arrows give brief descriptions of the approximations made, see references for further details. Note 
in particluar that the present formulation is a general one, and all existing formulations may be 
derived from it. 
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tensor in non-equilibrium thermodynamics |35j and thus is not negligible in general. 
Hence, for atomic liquids [2J, it has been assumed that it can be approximated by 
vdtp., where v is an arbitrary collision frequency. Although this resulted in a DDFT 
analogous to that previously derived for colloids in the high friction limit |39j . it 
is not clear that this is the correct approximation in general. For colloids with no 
hydrodynamic interactions, this term can be dealt with using a local-equilibrium 
approximation, or a Taylor expansion close to equilibrium, or considering the high 
friction limit [3]. However, the first two approaches are unsatisfactory for general 
systems which may not lie close to (local) equilibrium, whilst the high friction limit 
was not analysed rigorously. This high friction limit is the main objective of the 
present study. We will show that, in this limit, the term J dp (p <g> p — 1)/^ (r, p, r) 
is indeed negligible compared to p and the momentum distribution. 

It is worth noting that, if hydrodynamic interactions are neglected, the heuristic 
high-friction calculation made by Archer [3J produces the same DDFT as that derived 
by Marconi and Tarazona [3H1 SO] • We shall demonstrate that this still holds for the 
rigorous derivation. However, when hydrodynamic interactions are included, the two 
approaches do not lead to identical equations. Section [4] discusses in detail these 
differences. 

1.3. Toward a rigorous derivation of dynamic density functional theory. 

Our main result is that, for a system of N identical, spherically symmetric colloid 
particles, up to errors of C(e 2 ), where e ~ 7 _1 with 7 the friction constant for 
an infinitely dilute system (see Section [2| , the dynamics of the one-body position 
distribution p are given by 

M'.r)=-^V r -a(r,T), 



where a is the solution to a particular Fredholm integral equation (Theorem 4.1|, 
kg is Boltzmann's constant, T the absolute temperature and m the mass of the 
colloid particles. Furthermore, we show that if the one-body phase space distribution 
is written as a Hilbert [3T] (or Champman-Enskog |15j ) expansion /W(r,p,i) = 
/o(r, p, t) + e/i(r, p, t) + e 2 f 2 {v, p, t) + • • • then (in appropriate units) the first two 
terms are of the form p ( r j T ) e_ ' p ' ^ 2 and a(r, r)-p e~' p l I" 1 respectively. In particular, 
non-zero terms in the integral J dp(p ® p — l)/' 1 ' (r, p, i) are at most C(e 2 ). 

We note that the evolution equation takes the form of a continuity equation. In 
the framework of standard fluid dynamics one would expect a(r, t) — p(r,t)v(r,t), 
where p is the fluid density and v is the velocity field. Using the stan dard defini tion 



[3] v(r,t) : = g- 1 (r,t)/dpp/W(r,p,f), the Hilbert expansion d3.lt, Corollary 



3.3 



and Lemma 3.13 show that, up to errors of 0(e ), this interpretation holds. 

We now discuss the novelty of our approach and results. In previous work, the 
hydrodynamic interactions have been ignored. In this case, the leading-order term 



in the expansion in the inverse of the friction constant becomes linear (cf. (3.2a) 
where it is nonlinear), making the analysis significantly easier. As noted above, the 
full iV-body equations are also linear. In this work, we will consider two-body hydro- 
dynamic interactions, along with a two-body inter-particle potential, which require 
the approximation of the two-body distribution. As will be seen, a standard ap- 
proximation then leads to quadratic nonlinearities in the one-body equation, formally 
analogous to the quadratic nonlinearity of the collision operator in the Boltzmann 
equation [T7J [TTJ [S3]. For more general interactions, the nonlinearities will be of 
higher order; at least a priori, n-body interactions require n-body distributions. In 
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Section [4] we show that this heuristic argument does not actually hold when starting 
from the iV-body Smoluchowski equation; higher distributions are required. 

Let us now contrast the rigorous derivation of hydrodynamics from the Boltzmann 
equation (in the limit of small mean free path, or high collision frequency) - see 
e.g. the comprehensive review by Lebowitz et al [22] ■ The approach used therein, 
where time (and possibly space) are suitably rescaled and a Hilbert expansion is used 
to derive an infinite hierarchy of equations, which may then be solved to arbitrary 
order, is very similar in spirit to ours. The collision term is replaced by a term 
involving the hydrodynamic interactions, which has been much less widely studied 
than the Boltzmann collision operator. Determining the leading order term in the 



Hilbert expansion requires finding the null space of the collision term (see (3.2al) 



The full friction operator is a complicated integral operator, and determining its null 



space is non-trivial (see Lemma 3.2). In contrast, for the Boltzmann collision and 
self-friction operators it is straightforward to show that the null space contains only 
Maxwellians. Furthermore, in our situation, there are additional nonlinear terms due 
to the inter-particle potentials. However, due to these terms being independent of p, 
the momentum variable, and occurring with a higher power of the small parameter, 
they do not hinder the analysis in the same way as the hydrodynamic interaction 
terms. In addition, these nonlinearities affect our ability to control the evolution of 
the parts of the higher order corrections which lie in the null space of the operator we 
need to invert. Sections 3^ and [5] highlight in detail these difficulties. 



The structure of the paper is summarised as follows. In Section[2]we give a descrip- 
tion of the model, in both the original and rescaled timescales, state our assumptions 
and give an overview of the main result. In Section [3] we develop the solvability con- 
dition for the Hilbert expansion of the one-body distribution /W, which forms the 
basis for the proof of the main result stated stated in Section [IJ where we also discuss 
its relationship with existing formulations of the one-body Smoluchowski equation. 
In Section [5] we discuss the impact of our main result, including its application to 
the derivation of DDFT, and also describe a number of associated open problems. 
Appendix |A"[ contains proofs of the more technical lemmas of Section [3j 

2. Description of the model and statement of main results. We begin 
by considering the full equations of motion, in both position and momentum for a 
large number N of spherically symmetric colloid particles of mass to in a bath of a 
much larger number of much lighter particles. The interaction between the colloidal 
particles and the bath is modelled on the level of stochastic noise and the interaction 
between colloidal particles mediated by the bath is modelled by friction terms. The 
magnitude of these two effects is correlated due to a generalised fluctuation-dissipation 
theorem HQ HQ ES HI] . 

The evolution equations are 

N N 

3 = 1 j=i 

where, for r N = r u r 2 , ...,r N , r y G R 3x3 (r N ) and T = (r y ) € R^nxsn^n^ ig the 
friction tensor, which is positive definite, and in particular has a square root. 7 is the 
friction constant for a single isolated particle, and we are interested in the regime where 
7^1. The Wj(t) = (Wj(t),Wj(t),Wj(t)) T are mean zero, uncorrelated stochastic 
white noise terms and satisfy (it?™(t)) = and (w™(t)w™(t')) = Sjk6 nm S(t — t'). T 
is related to the strength of the stochastic white noise terms Wj via a generalized 
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fluctuation-dissipation theorem, namely A = VT. represents the force on particle 
i exerted by an external field and interactions with the other colloid particles. T is 
the absolute temperature, and k B is Boltzmann's constant. 

The motivation for this analysis is that in the high friction (overdamped, large 
7) limit, the momenta should reach equilibrium on a much faster timescale than the 
positions. In particular, wc are interested in times of order 0(7 _1 ) and so begin by 
rescaling the time variable as t = k B T / (m/yW (where t and r are the new and old times 
respectivel y), set X j = - V r V and define X 4 := Xi/(k B T) = - V r V = -V r U/k B T 
and e = Jk B T. '/mj -1 . This rescaling leads to the following system of equations: 

1 

r, = ~Pi 
e 

N N 

The constant in the time rescaling corresponds physically to Dq, the diffusion 
constant for a single isolated particle. The rescaling of X corresponds to measur- 
ing potential energy in units of the temperature. For e, we note that \/k B T /m is 
the average thermal equilibrium speed of a particle at temperature T, whilst 7 _1 is 
approximately the time required for the velocity distribution of the colloids to equi- 
librate. Hence e has units of length, and in order to produce a non-dimensional 
constant, it would be necessary to introduce another length scale. Such a scale is 
highly problem-dependent, and could for example be the typical length over which 
the external potential varies, the length of a finite box in which the particles are 
contained, a typical separation of colloid particles, or the size of the colloids. As 
such, we retain the dimensional parameter e and remark that the existence of a small 
parameter for applications must be checked on a case-by-case basis. 

The corresponding Fokker-Planck equation for the TV-body distribution function 

is 

N 



9 t /W(r",pV) + -f>- V ri /W(r^,pV) 



1 N r 1 

= j> E V p. ' [r y (r Ar )(p,+V Pj )/W(r w ,p JV ,i)]. (2.1) 



Here we have used the notation r" = (ri, . . . ,r„) and the analogue for p™. We also 
find it convenient to write dr^"" = dr n+ i . . .drjv and the analogue for p. In the 
above, f( N '(r N ,p N ,t) is the probability of finding each particle i at position r, with 
momentum p,j at time t. We note that the above equation is precisely the ./V-body 
Kramers equation, see e.g. [T§1 Hi] . 

It is clear that encodes a huge amount of information, and as such is very 
computationally demanding. We are not interested in the distributions of the positions 
and momenta of all the identical particles, but in the distribution of the average values 
of these quantities. For this reason we introduce the reduced probability distributions 

/W(r™,pV) := ^ ! J dr N - n dp N - n f^Hr N ,p n ,t). (2.2) 



Multiplying (2.1) by N and integrating over dr dp , all terms in the sums 



with i ^ 1 vanish and the evolution of the one-body distribution is given by 

N 

= ^ E / dr w - 1 dp^ 1 r 1 ,(r w )(p J +V Pj )/W(r w ,p^,i). (2.3) 

The difficulty in solving this equation lies primarily in the fact that the last 
two terms still involve f^ N ', the full iV-body distribution. In order to remove this 
dependence, and obtain a closed equation, it is necessary to make some assump- 
tions. Firstly we assume that the potential and friction tensor contain at most 
two-body interactions. We will show that this is equivalent to requiring knowl- 
edge of only f( 2 \ From previous studies on the derivation of DDFT (e.g. [14j), 
it is known that the full TV-body distribution function f( N ^ can be written as a 
functional of the one-body spatial distribution, and therefore so can f^ 2 \ We make 
the Enskog approximation to the two-body distribution, in particular assuming that 
/ (2) (ri, pi, r 2 , p 2 , t) = /W (n, pi., £)/W (r 2 , p 2 , t)g(T\,T 2 ), where the pair-distribution 
function g is assumed to be inde pendent of p and e. 



With this assumption, (^2| shows that / drdp/ (1) (r, p, t)g(r, r') = N - 1. The 
role of g is to describe positional correlations of the particles, such as finite-size ex- 
clusion effects. It is an intermediate level of approximation between the mean field 
approximation (g = 1) and the full 2-body distribution function f^ 2 '. We note that if 

= Sijl, the method outlined below allows a Smoluchowski equation to be derived 
even if g depends on pi and p 2 . It seems unlikely than an analogous result holds in 
general, as non-trivial momentum correlations on the two-particle level would prevent 
the system from thcrmalising to a Maxwellian momentum distribution. 

The assumption that g is independent of e is known not to be valid in general; we 
have only that g = g(r, r'; [p]) (g is a functional of p), where p(r, t) = J dp/W (r, p, t). 
This assumption does hold if the only contribution to p comes from the zeroth order 
term in /W. Even so, as we will discuss in Section pi if we expand g in a power 
series in e, g = go + egi + . . . , then the derivation changes only at the e° level, and 
we recover an analogous Smoluchowski equation. Such an approximation is standard 
in the physics literature and, as suggested by the name, was first proposed by Enskog 
[lOj and revised by Van Beijeren and Ernst |70j to ensure consistency with irreversible 
thermodynamics. See also [53]. 

To summarize: 
• Assumption 1 [Pairwise additive potential] 

N 



V(r N ,t)=J2Vi(r h t) + lj2 V ^ r ^ 



2 

i=i i^j 
• Assumption 2 [Pairwise additive friction] 

r y = Sijl + f y = <%(1 + Z i( r *^)) + (! - %)Z 2 (ri,rj) 

= % J2 (iv^T 1 + Z i( r *> r ^)) + i 1 ~ %)Z 2 (r,,r,), 
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with the Zij symmetric 3x3 matrices. 
Assumption 3 [Enskog approximation] 

/ (2) (ri,pi,r2,p a ,t) = /«(r l! p 1 ,t)/«(r 2 ,p 2 ,t) 9 (r 1 ,r 2 ). (2.4) 



Making these assumptions, we now calculate the two remaining terms in (2.3 1. 
Using standard symmetry arguments gives: 

N J dr N ~ 1 dp N ' 1 V ri V(r N )-V Pl f^ N '>(r N ,p N ,t) 

= V ri V 1 (r 1 ,t)+ f dr 2 dp 2ff (r 1 ,r 2 )/( 1 )(r 2 ,p 2 ,i)V ri U 2 (r 1 ,r 2 )l ■ V Pl f^ (r llPu t). 



and 



= V pi • [(Pl + Vpj/W^p!,*) 

+ J dr 2 dp 2ff (r 1 ,r 2 )/( 1 )(r 2 ,p 2 ,i)Z 1 (r 1 ,r 2 ) x (pi + VpJ/^n.Pi,*) 
+ J dr 2 dp 2ff (r 1 ,r 2 )Z 2 (r 1 ,r 2 )(p 2 + V P2 )/( 1 )(r 2 ,p 2 ,i) x / (1) (ri,Pi,t) 



Hence, we have: 

Proposition 2.1. Under Assumptions 1, 2 and 3, the evolution of the one-body 
reduced distribution satisfies 



5 t /« (r, p, t) = - [ - p • V r + V r Vi(r) • V p 

+ J dr'dp'f( 1 Hr',p',t)g(r ) r')V T V 2 (r ) r').V p \f^(r,P,t) 
[(P + V p )/«(r,p,t) 

+ f di/dp'sCryjZiCryj/W^p',*) x (p + V p )/«(r,p,i) 
+ 1 dr'dp'. 9 (r,r')Z 2 (r,r')(p' + V p O/ (1) (r',p',i)x/«(r,p,i) 
- + +Ar (/W,/«)], (2.5) 



where 



Cof = V p -(p + V p )/(r,p,t) 
£i/=[-p-V r + V r F 1 (r)-V p ]/(r,p,t) 



(2.6a) 
(2.6b) 

dr'dp'.g(r, OZ^r, r')/(r', p', i) x (p + V p )/(r, p, t) (2.6c) 
+ V p .y dr'dp'. 9 (r,r')Z 2 (r,r')(p' + V P 0/(r',p',t)x/(r,p,i) 
dr'd P 7(r',p',t)g(r,r')V r y 2 (r,r') • V p /(r,p,t), (2.6d) 



and for ease of notation we have omitted the (r, p, t) dependence of the functions on 
the left hand sides of (2.6a)-(2.6d). 

It is noteworthy that, although they are quadratic, the non- linear terms are not 
symmetric in the two arguments. The first argument has been taken to be that inside 
the integral. 

In the following we will assume that f^ is bounded and positive^] and that all 
functions are sufficiently regular and have sufficient decay at infinity for operators and 
integrals to be defined. 

To state our main result, we recall that the position distribution, which is the 
object of interest in the single-particle Smoluchowski regime, is defined by p(r, t) = 
J dp/W (r, p, t). We will show that its evolution equation is given by 

Theorem 2.2 (Smoluchowski equation). Under suitable assumptions on f^ 1 ', 
Uj , Zij, j = 1,2 (see Theorem 4-1 K U P t° errors of 0(e 2 ) the dynamics of the one- 
body position distribution are given ( in the original timescale ) by: 

M^) = -^V r .a(r,T), 

where a(r, r) is the solution to 

a(r,r) + J dr' 5 (r, r')p(r', r )Zi(r, r') x a(r, r) + p(r, r) J dr' 5 (r, r')Z 2 (r, r')a(r', r) 

= - V r + ^(v r C/i(r,r)+ /drV(r',r) 5 (r,r')V r [/ 2 (r,r'))]p(r,r). (2.7) 



3. The Hilbert Expansion. We now expand in powers of e as: 

oo 

/«(r,p,i) = ^e n / n (r,p,i). (3.1) 

Due to the singular nature of the problem, we do not expect such a regular perturba- 
tion expansion to converge uniformly. The expansion should be valid only for times 
t 3> e, and not for shorter times, i.e. we expect there to be a boundary layer in time 
of size 0(e). Since we are interested in times much larger than e, interest lies in the 
leading order terms; one would then hope to be able to truncate the series and prove 

is non-negative by definition, but may be zero, e.g. if there are excluded areas of the phase 
space due to confining potentials. It is clear that it may be made positive with arbitrarily small 
errors. 
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suitable bounds on the remainder term, as in [27]. We also assume that such an 
expansion then converges, in particular that the /„ are sufficiently well-behaved in r 
and p. 



Inserting (3.1 1 into the evolution equation (2.5 1 and collecting powers of e gives 



the hierarchy of equations 



£o/o+.A/o(/q,/o)=0 (3.2a) 

4>/i +CJ +A/- (/ , A) +AA (A,/ ) +M(/o,/o) = (3.2b) 

4)/2+A/l +M)(/2,/o) +M>(/0,/ 2 ) +AA (/i,/i) 

+M(/i,/o) +M(/o,/i) = fit/o (3.2c) 



i+j— n— 1 



We now solve these equations order- by-order. First, to solve (3.2a I, we need to 



determine the null space of £q ■ +Afo(-, •)■ However, before we do so, the following 
lemma will be useful. It is essentially a result of the positive-definiteness of T: 

Lemma 3.1. For v(r, p, t) an arbitrary vector such that the integrals below exist, 



and f satisfying (2.4-), there exists 8 > such that 



J drdpdr'dp7(r, p, t)f(r', p', %(r, r') 

x [v(r,p,t) ■ (jv^l + Z 1 (r,r'))v(r,p,t) +v(r,p,i) • Z 2 (r, r')v(r', p', t)] 
><5 /" drdp/(r,p,i)|v(r,p,i)| 2 . 



JVi particular, the result holds when f is chosen to be either fo or f^K 
Proof. See Appendix | A. 1 1 □ 



solution to (3.2a) 



We now proceed with the analysis of ( 3.2a |— ( 3.2c), beginning by determining the 



3.1. Solution of the e 2 equation. In this section we find the solution fo of 



(3.2a|: 



Lemma 3.2. For f satisfying {2.4), the null space of Cof +-A/q(/, /) consists of 
functions of the form f(r,p,t) = exp ( — |p| 2 /2) 0(r, t). 



Proof. We begin by assuming that / satisfies (2.4) and Cof + A/"o(/, /) = 0. We 
define <j) by /(r, p, i) =: e~l p l"/ 2 (f>(r, p, t) and note that is positive. We therefore 
have m<£[£ /+M>(/,/)] = 0. Note that /(r,p,t) = dr'dp'/( 2 )(r, p, r', p', t) 

and so 



^o/(r,P,t) = 7v L T V P' J dr'dp'«?(r,r')/(r',p',i)(p + V p )/(r,p,i), 
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Using this reformulation along with the definition (2.6c), the notation Zi(r,r') := 
ty3t1 + Zi(r, r') and integrating over r and p gives 



0= / drdpln<£(r,p,t)V p - J dr'dp' ff (r, r') [/(r', p')Z x (r, r')(p + V p )/(r, p, t) 
+ /(r, p, t)Z a (r, r')(p' + V p ,)/(r', p', t) 



drdpdr'dp'e- |p| 12 e~ lp ' 1 /2 g(r,r') 



^ V p ^(r,p,t) 
(r,P,*) 



(r', p', t)Z x (r, r')V p( J(r, p, t) + 0(r, p, t)Z 2 (r, r') V p ^(r', p', f) 



where we have used integration by parts and Fubini's theorem, along with the identity 
(p + V p )[e-lpl 2 / 2 0(r,p,i)] = c-IpI 2 / 2 V p 4>(r,p,t). 
Letting v(r, p, t) := V p <^>(r, p, t) /<j)(r, p,t) we have 







drdpdr'dp'e Hp|/2 e^ 1 * 5 ' 1 /2 <p(r, p, i)0(r', p', i)s(r, r') 

v(r, p, i) • Zi(r, r')v(r, p, i) + v(r, p, i) • Z 2 (r, r')v(r', p',t) . 

which shows 



3.1 



Since /(r, p, t) = exp (— |p| 2 /2)</>(r, p, t) > 0, we may apply Lemma 
that 

= J drdp/(r,p,i)|v(r,p,i)| 2 = / drdp ' 2 <\T x (r, p, t)\ V p 0(r, p, t)\ 2 . 

Hence, since 4> is bounded, the integrand is zero if and only if V p 0(r, p, t) = and 

the result holds. □ 

This result gives an explicit form for the p-dependence of /o- 

Corollary 3.3. The zeroth order term in the e-expansion of f^ is given by 



/o(r,p,t) 



1 



(2tt) 3 / 2 



e = Po (r,t) =: Z L c * p (r,t). 



Proof. Follows immediately from (3.2a) and Lemma 3.2 □ 

3.2. Solution of the e _1 equation. In order to find f\ from (3.2b I, we rewrite 



it as 



4>/i +AA (/o,/i) + JV (/i,/o) = -A/o -M(/o,/o). 



(3.3) 



The first point to note is that for known /o, this is a linear-operator equation for /i. 
Although this can be seen in an abstract sense from the non-linearities being of a 
quadratic nature, we will require the explicit form of the operator: 

Lemma 3.4. For No as in (2.6c), /o(r, p, t) as given by Corollary 3.3 and 
arbitrary /(r, p, t), 



Wo, /) = V p • [ / dr'dp'«?(r, r')/ (r', p', ^(r, r') x (p + V p )/(r, p, t) 
A/o(/,/o) = -/o(r,p,t)p- J dr'dp' ff (r,r')Z 2 (r,r')(p' + V P 0/(r',p',t). 
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Proof. For A/o(/o, /) note that (p + V p ) e IpI 2 / 2 = and hence the second line in 



(2.6c) gives zero and the first result follows. For Ao(/, /o), the same argument shows 



that the first line in (2.6c) is zero. For the remaining term, the result of the integral is 
a vector which depends only on r, and for any such vector z, V p - [z(r)/(r, p, i)] = z(r)- 
V p /(r, p, t). The result then follows from the identity V p /o(r, p, t) = — p/o(r, p, t). 
□ 

COROLLARY 3.5. For known fo, and arbitrary f, Cf :— Cof + A/o(/o,/) + 
A/"o(/, fo) is a linear operator on f. 

As is customary for problems of this form, it is convenient to work in the L? space 
weighted by the inverse of the invariant measure of C In this case, the inner product 
is defined as 

(/(r, P, t), /(r, p, t)) f -i := J drdp/^O, p, i)/(r, p, t)f(r, p, t). 

We denote this space by L 2 i. In this weighted space, C is self-adjoint and has 

to 

compact resolvent, allowing us to apply Frcdholm's theory: 
Lemma 3.6. C is self-adjoint in L 2 i. 

Proof. See Appendix |A.2| D 

Lemma 3.7. The resolvent of C is compact in L 2 1 . 

to 

Proof. See Appendix |A.3| □ 



Since (3.3) may be rewritten as Lf\ = — £i/o — A/i(/o,/o), to determine the 
solvability condition we must determine the null space of C* . 

Lemma 3.8. The null space of C (and of C* ) contains only functions of the form 
/(r,p,t)=cxp(-|p| 2 /2)0(r,O. ' 

Proof. The proof is analogous to that of Lemma [3~2"1 but is more straightforward 



as the operator is now linear and we are already working in the appropriate weighted 
space. As in the proof of Lemma 3.2 we have 

^o/(r,p,<) = 1 v L T V p - j dr'dp'.g(r,r')/o(r',p',i)(p + V p )/(r,p,t), 

and hence, setting Zi(r,r') := t^ttI + Zi(r, r') we have 

£/(r,p,i)=V p - j dr'dp' ff (r,r')/o(r',p',i)Z 1 (r,r')(p + V p )/(r,p,i) 

+ g(r, r')/o(r, p, t)Z 2 (r, r')(p' + V P 0/(r', p', t) = 0. 
Taking the L 2 , inner product with / gives 

to 







J drdp/o" 1 (r,p,t)/(r,p,t) 

x V p .y dr'dp'[. 9 (r,r')/o(r',p',t)Z 1 (r,r')( P + V p )/(r,p,0 

+ g(r, r')/o(r, p, t)Z 2 (r, r')(p' + V P 0/(r', p', t) 



Integrating by parts, using the identity V p / x (r, p, t) = p/ 1 (r,p,t) and Fubini's 
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theorem gives 

= J drdpdr'dp7 - 1 (r,p,i)(p + V p )/(r,p,t) 

• g(v, r')/o(r', p', t)Zi(r, r')(p + V p )/(r, p, *) 

+ g(r, r')/o(r, p, t)Z 2 (r, r')(p' + V p ,)/(r', p', t) 

Setting v(r,p,t) := /q^Oj p, t)(p + V p )/(r, p, f) then gives 

= ^ drdpdr'dp7 (r,p^)/o(r',p',t) ff (r,r')v(r,p,t) 

• [Z 1 (r,r')v(r,p^) + Z 2 (r,r')v(r',p')]. 

Hence, by Lemma |XT| we have 

0- J drdp/ (r,p,i)|v(r,p,t)| 2 = | drdp/^r, p, t)|(p + V p )/(r, p, t)\ 2 . 

Since / is bounded and positive, we must have (p + V p )/(r, p, i) = 0. Writing 
/(r,p,t) = Z~ x exp(— |p| 2 /2)</>(r, p, t) then gives V p </>(r, p,f) = and the result 
follows. □ 

Determining the explicit solvability condition finally requires the explicit calcula- 



tion of the right hand side of (3.3) 



Lemma 3.9. For Ci and A/i as in (2.6b) and (2.6d), and fo as in Corollary 3.3 



&fo = -Z^c-^-p- [V r + V P Vi(r)]po(r,i) 



and 



M(/o,/o)(r,p,t) - -Z-'e-^ Po (r,t)p- / dv'p (r', t)g(r, r') V r V 2 {v, r') 



Proof. This is a simple calculation: For £i we have 
£i/o(r, p, i) = [ - p • V r + V r ^(r) • V p ] \z~ x po(r, t) 

= Z- 1 e-^ i [-p. V r po(r,t)] +V r ^(r). [ - p] Z" 1 e~^ Po (r, t), 

and the result follows. 

The result for M\ follows from the two identities J dp/o(r, p,t) = po(r,t) and 
V p / (r,p,<) = -pc-|pl 2 /2 po ( r;t ). □ 



Recall that we are trying to solve (3.3), and have shown that in L 2 i, C is self 

I I Jo 

adjoint with compact resolvent, and has null space elements exp(— |p| 2 /2)0(r, t). In 



order for (3.3) to be soluble, we therefore require, by the Fredholm alternative, that 
its inner product (in L 2 ,) with any element of the null space of C* is zero. 

Note that (e~' p ' I" 2 <j)(r, t), = (po( r , t)^(r, t), f) and we therefore require 

that the integral with respect to p of the right hand side of (3.3) is zero. This is 
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an easy corollary of Lemma |3.9| since the p dependence of both terms is of the form 
exp(— |p| 2 /2)p, which integrates to zero. 

Corollary 3.10. —Cifo — A/i(/o,/o) is orthogonal to the null space of C* and 



thus (3.3) always has a solution. 



Since we now know that (3.3) is soluble, we can invert its left hand side. The 



standard approach would be to expand in a basis of the eigenfunctions of C. However, 
since 2*1 and Z 2 (which enter C through Ao) are unknown, we expand in a basis of 
products of generalized Hermite polynomials multiplied by a Maxwellian, which are 
eigenfunctions for the case Zi = Z 2 = 0. This turns out to be sufficient as we do not 
need to explicitly invert C 

Definition 3.11. We define the basis of L 2 (K 3 , e~l p l / 2 ), for p = (pi,P2,P3) T , 

fn,a(p) := H ai ( Pl )H a2 (p 2 )H a3 (p 3 ), 

where n G N. a = (ai,a2,o.3) T , Oj € N, |a| := ^ = n, and H n are the standard 
one- dimensional Hermite polynomials. 

Since M 3 is a product space and the Hermite polynomials form an orthogonal basis 
of L 2 (R, e~P 2 / 2 ), it is clear that the F„, a form an orthogonal basis of L 2 (M 3 , e _|p|2/2 ) 
(see e.g. [53]). We now show that C preserves n, the degree of the Hermite polynomial, 
when applied to e~' p ' I" 1 -P„. a (p)- Hence, with a slight abuse of notation, we find 

£-l [c -|p| 2 /2 p n a{p)] e Span {{ e -| P | 2 /2 p n h{p) | | b | = n} [J e -|plV2 } ; 

where the span runs over coefficients in r and t, and we note that e~l p l I" 1 is the 
p-dependent part of the kernel of C 

Lemma 3.12. Denote the expansions of Cf and f by 

oo 

Cf=: e-l p l 2 / 2 ]T Y, 7n,a(r,t)F„, a (p) 

n— |a|— n 

oo 

/=: e -l p l 2 / 2 £ ]T 7 n,a(r > *)P„,a(p). 

n=0 |a|=n 

Then, for each fixed n > 1, 7n, a (r,i) = for all a if and only if "f n ^ a (r,t) = for all 



Proof. See Appendix |A.4| D 

We are now in the position to determine the explicit form of f\ and thus solve 



(3.2b I 



Lemma 3.13. /i(r, p,t) = [a(r, t) ■ p + i>(r, t)]Z 1 exp(-|p| 2 /2), where a(r,t) is 
given by the solution of 

a(r,t) + J dr' ff (r,r')p (r',0Zi(r,r') x a(r,t) + p (r,t) / dr' 3 (r, r')Z 2 (r, r')a(r', t) 
= - Vr + V r Fi(r,t)+ [ dr'p (r',t)g(r,r')V r V 2 (r,r')]p (r,t). 



Proof. By (3.3 1 and Lemma 3.9 we have 



£/i = -£/o-M(/o,/o) 

= V r + V r V 1 (r,t) + J dr , po(r , ,%(r,r , )V r F 2 (r,r , )]/9o(r,f)-p^ 1 e-l p l 2 / 2 

=: a(r,<) •pZ- 1 e- |p|2/2 . (3.4) 
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Hence, by Lemma 3.12 and the definitions P\, ej = Pj for &j the standard unit vectors, 
it follows that /i(rTpTT) = [a(r, t)-p+^(r, t^Z' 1 exp(- |p| 2 /2) for some a. Evaluating 
each of the terms of Cf\ then gives, firstly from (2.6a I, 

A>/i = V p • [(p + V p )[a(r, t) ■ p + V(r, t)]Z~ x ' 2 ] 
= V p ■ [Z- 1 e-IPl 2 / 2 V p [a(r, t) ■ p + tf(r, t)]] 
= V p • [Z" 1 e-IPl 2 / 2 a(r, t)] = -Z' 1 e^' 2 p • a(r, t). 



Using the explicit form of Ao(/o, /i), as given by Lemma 3.4 gives 
Wo,/i) = V p - f / dr' 3 (r,r')p (r',t)Z 1 (r,r') 



x (p + V p ) [[a(r, t).p + V^(r, t)\Z~ x ^ I 2 } 
= V p • [ J dr' 3 (r, r')po(r', t)Z x {v, r') x a(r, ijZ" 1 ^ ' 2 
For a general matrix Z(r), we have 

3 

V p - [Z(r) aL (r,t)Z- 1 e~M 2 / 2 ] = £ 9 Pi [Zy(r, l^Z' 1 e'^^ 2 ] 

3 

= - £ ^(r.t^a^-^-IPl 2 / 2 =-p-Z(r,t)a(r,i)Z- 1 e-|Pl 2 / 2 , 

and hence 

A^o(/o,/i) = -^ 1 e- |p|2/2 p- J dr / <7(r,r / )p (r , ,t)Zi(r,r / )a(r,t). 

For the third term, we take A/"(/i,/o) as given by Lemma 3.4 and note that 
(p' + V P 0/i(r',p',i) = a(r',t)e-|p'l 2 / 2 , giving 



AA (/i, /„) = -Z" 1 dPl 2 / 2 p (r, i)p ■ / dr'dp' 5 (r, r')Z 2 (r, r')a(r', tjZ" 1 e"^ 
= -Z- 1 e-|Pl 2 / 2 p (r,<)p- 1 dr' 5 (r,r')Z 2 (r,r')a(r',t). 
Collecting the three terms gives 

a = -a(r, t) - J dr'. 9 (r, r> (r', t)Z x (r, r') x a(r, t) 
-po(r,t) J dr'. 9 (r,r')Z 2 (r,r')a(r',i), 
and the result follows by (3.4). □ 



3.3. Solution of the e equation. We have, by Corollary 3.3 and Lemma 3.13 



that 



/o(r,p,i) = Z- 1 exp(-|p| 2 /2)p (r,i), 

/i(r,p,t) = [a(r,t) • p + ^(r, t)]^ 1 exp(-|p| 2 /2). 
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We now show that the evolution equation for po is given by the solvability condition 
for the equation corresponding to e°, namely (3.2c). We begin by rewriting (3.2c) as 



th = A/i +AT (fiJi) +M(/i,/o) + M(/o, A) -dtf . 



(3.5) 



We then require that the right hand side is orthogonal to the null space of C* in L 2 1 

Jo 



which, by Lemma |3.8[ is equivalent to it being orthogonal to constants in p in the 
unweighted space L 2 . Hence, using the divergence theorem and the explicit forms of 
Mo and N\ as given by (2.6c) and (2.6d), the requirement reduces to J dp(£i/i — 
dtfo) = 0, and we need only calculate these two terms. However, for completeness 
and later use, we also calculate the A/q and J\f\ terms: 



A/" (/i,/i 



_ Z -l e -|p| 2 /2 p 
4- p"Ip| 2 /2 ' 



Lemma 3.14. Ford, No andNi as in (2.6b)-(2.6d), f as in Corollary 3.3 and 
fi as in Lemma \S '. 1 3\ we have 

dh = Z- 1 c-IpI 2 / 2 [ - p • V r a(r, i)p - p • V r VXr, t) 

+ V r V x {v) • [a(r, t) - p(p • a(r, t) + ^(r, t))] 

M(/o,/i) - Z^e-IPl 2 / 2 J dr'p (r',i) 5 (r,r')V r F 2 (r,r') 

• a(r,t)-p(p-a(r,t) + V(r,t)) 

M(/ 1 ,/o) = -^- 1 c-|Pl 2 / 2 p (r,i)p. j drV(r',<).g(r,r')V r y 2 (r,r') 

dr' 5 (r, r')^(r', f)Z x (r, r')j a(r, t) 
dr' ff (r,r')Z 2 (r,r')a(r',i) 

• a(r,t)-p(p-a(r,t)+V(r,<)) 

Proof. We begin with C\fx, which is given by 

£i/i(r, p, t) = [-p • V r + V r Vi(r) ■ V p ] [Z' 1 e^^ 2 (a(r, t) ■ p + 4>(r, t))] 
Simple calculations show that 

Vp^e-IPl 2 / 2 (p • a(r,i) + ^(r,t))] = Z^ 1 e^l 2 / 2 [-p + V p ](p • a(r,i) + ^(r,t)) 
= Z- 1 e-IPl 2 / 2 [a(r, i) - p(p • a(r, t) + ^(r, <))] (3.6) 
-p • VriZ- 1 e- |p|2/2 ^(r, i)) = -Z" 1 e- |p|2/2 p • V r V(r, t) 

This gives three of the terms. The calculation of the fourth is easier in coordinates. 
Using that p • V r = X)j=i P»^r< and p • a = 5Z J=1 Pjdj gives 

-p ■ V^Z- 1 e- |p|2/2 p ■ a(r, t)] = -Z" 1 e - |p|2/2 p • V r [p • a(r, *)] 

3 

= _ z -i e -\ P \ 2 /2 Pl d r A Pj a 3 (r,t)} 
= -Z^c-^^p- V,a(r,t)p, 
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which gives the result for C\j\. 

The expressions for A/i(/o,/i) and A/i(/i,/o) result from the trivial identities 
Jdp7o(r',p',t) = p (r,t), /dp7i(r',p',t) = ip(r,t), V p / (r,p,<) = -p/ (r,p,t) 
and (3.6). 

Finally, for A/"o(/i,/i) we use the trivial identity J dp' fi(r' 7 p', t) = ip(r,t), and 
that (p + V p )/i(r, p, t) = Z- 1 e"IPl 2 /2 Vp [p . a(r, t) + yj(r, t)] = Z~ x e"lpl 2 /2 a (r, t), 
and so / dp'(p' + V p /)/i(r, p', t) = a(r', t). Then 



M)(/i,/i) = V p 



dr' 5 (r,r')V>(r',t)Zi(r,r') x a(r, t)Z^ e~ |p|2/2 
dr' ff (r, r')Z 2 (r, r')a(r', t) x Z" 1 e~ |p|2/2 (p • a(r, t) + ip(r, t)) 



both terms of which are of the form V p • v(r, t)(j>(r, p, t), where v is a vector. Using 
V p .v(r,^(r,p,i) = v(r,t) ■ V p 0(r,p,t), V p e-I p l 2 / 2 = -pe^ 2 / 2 and (Si 
completes the proof. □ 

Recall that we wish to solve ( |3.5[ ), and require that J dp(£i/i — dtfo) = 0. The 
identities / dpZ^ 1 e~l p l 2 / 2 = 1, / dpZ" 1 e~ |p|2/2 p i p i = show that 

/ dpZ- 1 e-l p l 2 / 2 [a(r,<)-p(p.a)]=a(r,t)- / dpZ" 1 e^ 2 / 2 p^^a,(r, i) 

^ i=l 
= a(r,i) - a(r,t) = 

(which also follows from these terms resulting from the V r Vi • V p term in L\ and the 
divergence theorem) and 

J dpZ- 1 e-l p l 2 / 2 p-V r a(r,i)p = J dpZ' 1 e^ 2 / 2 £ tw^a^r, t) = V r -a(r, t), 

and hence J dp£i/i = — V r -a(r,t). Since f dp9 4 / (r, p, t) = d t po(r, t), the solvability 
condition becomes 



d t p (r,t) = -V r • a(r,i), 



(3.7) 



which is precisely the equation describing the one-body position distribution evolution 
for the Smoluchowski equation, as given in Theorem |4.1| 



3.4. Solution of the e 1 equation. We now demonstrate that ip(r, t) = if 
i/)(r,0) = 0. This should result from the solvability condition for the e 1 equation, 
which has the form 

-Ch = A/a +M)(/ 2 ,/i) +M)(/i,/ 2 ) +M(/ 2 ,/o) + M(/o,/a) + M(/i,/i) - fit/i- 

Since once again the Ao and A/i terms do not contribute to the Fredholm alternative 
calculation, we have 



/" dpOCi/a - ft/i) = 0. 



From (3.5 1, we have that 



/ 2 = (-£) - 1 [A A + M) (A , A ) + M (A , /o ) + M (/o , A ) - ft /o] , (3.8) 
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with dtfi = Z^ 1 c"' p ' I" 1 [p • <9 t a(r, t) + dtip{r,t)] and the remaining terms given by 
Lemma 13.141 

For the £1/2 term, we have £1/2 = [— P • V r + V r Vi(r) • V p ]/2(r, p, t), and, by 
the divergence theorem, the second term vanishes upon integration. Hence we are 
interested only in 

J dp [Z' 1 e-IPl 2 / 2 d t t/j(r, t) + p ■ V r / 2 (r, p, *)] = d t i/>(r, t) + J dp p ■ V r / 2 (r, p, t). 

Since H\ ej (p) = pj, the only terms from / 2 which contribute to the integral are of the 
form p-a 2 (r,i)iT _1 e~' p ' 2 / 2 , i.e p-Vi(f2) where V\ is the projection onto p, i.e. V\f = 
J dpp/(r,p,t). By (3.8) and Lemma 3.12 it therefore suffices to consider only terms 
of the form -pa 2 (r, t)Z~ x e -|pl 2 /2 mA/i+ATo^, /o)+M(/o, fi)-d t fo, 

i.e. a 2 (r,t) = -V^Cxh +M(/i,/o) +M(/o,/i) - flt/o). 

By Lemma 3.14 and d t f (r,p,t) = Z- 1 c^ 2 / 2 d t p {r 7 t), we have 

Vx (-£/ 2 ) = Pi ( dh + Mo {h , h) + Mi {h , /o) + M (/o, /1) - dtfo) =: -a 2 (r, t) , 
with 



a 2 (r,i) := [V r ^»(r, i) + ^(r, t)V r 7 x (r) + / dr'p (r, %(r, r')V r V 2 (r, r')^(r, i) 
+ p (r,t)( / drV(r',t).g(r,r')V r y 2 (r,r')) 
J drV(r',t) 9 (r,r / )Zi(r,r / ))a(r,i) 

dr'« ? (r,r')Z 2 (r,r')a(r',i))^(r,Ol, 



where a(r, t) is given by Lemma 3.13 Thus, by (the proof of) Lemma 3.13 

7>i/ 2 (r,p,i) =a 2 (r,t) 

with a 2 the solution of 

-a 2 (r,i) = a 2 (r,t) + / dr'p(r, r> (r', i)Z x (r, r') x a 2 (r, t) 
+ p (r,t) J dr' ff (r,r')Z 2 (r,r')a 2 (r',i) 
Hence, for the e 1 equation to be solvable, 
= d t V>(r,i)+ /" dp p- V r Pi/ 2 (r,p,t) 
= 9^(r, i) + J dpZ- 1 e-IPl 2 / 2 p V r [p a 2 (r, t)) = d t ^(r, t) + V r ■ a 2 (r, f) 



or 



dtVM) = -V r -a 2 (r,t). 



(3.9) 



To ensure that ip(r, t) = 0, we first note that ^(r, 0) = is equivalent to assuming 
that the initial condition f^'(r, p, 0) is independent of e. For this to hold for all t, it 



19 



is necessary to show that (3.9| is dissipative (or, since dt J dr0(r,i) = 0, that (3.9 1 is 
non- negativity preserving) . 

The proof in the linear (V2, Zi, Z2 all zero) case is trivial, as it turns out that if> 
and p satisfy the same equation. Thus, since the Smoluchowski equation for p must be 
non- negativity preserving, so must the equation for ip. The proof in the general case 
is complicated both by the equations for p and ip not being identical (due to the non- 
linear terms) and by needing to prove dissipativity results for the resulting non-linear 
operators. In general (for Z 2 7^ 0), the equations are not even explicit as one needs 
to solve the Fredholm integral equations for a and a 2 . However, since the full friction 



tensor is positive-definite, one would expect (3.9| to be a parabolic PDE and so, for 



potentials and hydrodynamic interaction terms with sufficient bounded derivatives, 
the result should follow from standard PDE theory, see e.g. [75]. We therefore assume 
that Vi,V 2 , Zi and Z 2 are such that that if ^(r, 0) = then ip(r, t) = for all t > 0. 

4. The Smoluchowski equation. We are now in a position to state our rigor- 
ously derived Smoluchowski equation. For ease of comparison with existing results, 
we return to the original scalings of time and potentials. 

Theorem 4.1 (Smoluchowski equation). Suppose /W(r,p, 0) = / (r, p, 0) = 
^-ig-|p| / 2 p o (r,0) is independent of e. Suppose further that po(r,0), Uj and Zj, 
j = 1 , 2 a re such that the solutions of (2.5), [3.1) and (3.9) exist for times [0, to] and 
that (3.9) is non-negativity preserving. Then, up to errors of 0(e 2 ), the dynamics 
of the one-body position distribution are given (in the original timescale) for r G 
[0,mj/(k B T)t o } by 



d T p(r, t) 



_ k B T 
7717 



V r -a(r,r) 



where a(r, r) is the solution to 



a(r,r) + J dr' 5 (r, r')p(r', r)Z x (r, r') x a(r, r) + p(r, r) J dr'«?(r, r')Z 2 (r, r')a(r', r) 
= - V r + ^(v r C/ 1 (r,r)+ f dr'p(r',r). g (r,r')V r (7 2 (r,r'))]p(r,r). 



(4.1) 



Proof. The evolution equation is given by (3.7), which is the solvability condition 



for (3.2c), and a(r, t) is given by Lemma 3.13 Returning to the original timescale 



introduces the factor of p, = k B T / [m^f) in the right hand side. We also repl ace V by 
its original value of [//(fc^T) where X 4 = — V ri U(r N ). The conditions on (3.9) and 
the initial condition ensure that, using the notation of Lemma 3.13 i^(r, t) = for all 
times. Hence p(r, r) = po( r , T ) + ^(e 2 ). □ 

We note here that the assumptions on the initial condition and on the existence 
of solutions are analogous to those made for the Boltzmann equation, see e.g. [27] . 
We expect that proving such assumptions hold for a physically interesting range of 
potentials and friction tensors would be a formidable problem in its own right, and is 
beyond the scope of the present study. An analysis of the corresponding problem for 
the Boltzmann equation is given in |22j . 

To demonstrate the connection to existing formulations, we assume that Z 2 = 0, 
which allows us to find a explicitly. We then have: 

Corollary 4.2. 



Under the same assumptions as in Theorem 4-1, «/Z 2 
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0, the 



one-body position dynamics are, up to errors of 0(e 2 ), governed by 



d T p{v, t) = V r • (D(r, r) [V r p(r, r) + ^(r, r) V r ^(r, r) 

+ drV 2 )(r,r',r)V r F 2 (r,r')]), (4.2) 

where we have defined p^ (r, r', r) := p(r, r)p(r', r)g(r, r', r), as it would be for the 
Enskog approximation, and the 3x3 diffusion tensor D is given by 

D(r,r) = — [l+ / dr' 5 (r,rXr',T)Z 1 (r,r')j _1 
m7 L y J 



Corollary |4.2| gives a one-body Smoluchowski equation with a novel form for the 
diffusion tensor D. As is clear from the notation, D(r, r) not only depends on the 
position but also on the time. This time-dependence is present through the time- 
dependence of p, against which the two-body terms must be averaged. 

One obvious question is whether D is positive definite. A simple calculation shows 
that this is indeed the case. Note that 1 + ^i( r i; r i) i s positive definite (since it 

is a principal minor of T, which is positive definite). Hence, for any v(r, r), we have, 
for some S > 0, 



v(n, t) + z i( r i> r i)]v(ri,r) > 5|v(n,T) 



j¥1 



v(ri,r) • J dr 2/ 5(r 2 ,T) 3 (ri,r 2 )[ 1 Y ;L T l + Zi(n, r 2 )]v(r 1( r) > <5|v(ri, r)| : 



where the proof is virtually identical to that of Lemma pO| except we do not integrate 
over 17. Since J dr 2 p(r 2 , r)g(r, r 2 ) = N — 1, and &bT, m and 7 are positive, this is 
equivalent to D _1 , and hence D, being positive definite. 

We now compare our result with that derived by Rex and Lowen [SH (5)-(8)]. As 
demonstrated in Figure |L~H their Smoluchowski equation is derived from the ./V-body 
Smoluchowski equation for pairwise additivity of both the potential (our Assumption 
1) and diffusion tensor. The second assumption is analogous to our Assumption 2, 
but not equivalent, as the inverse of a matrix (recall TD = fc^T/ml) with pairwise 
terms need not contain only pairwise terms. However, there are situations where the 
two assumptions are essentially equivalent, such as in a diffuse colloid system. The 
underlying assumption then is that there exists an additional small parameter, say A, 
with 1 > A > e and such that Zi = 0(A). Then, up to errors of C(A 2 ), D(r,r) = 
— / dr'g(r, r')p(r', r)Zi(r, r')]. and thus the diffusion tensor is a two-body one. 
We note that the analogue of Assumption 3 is p^ (ri, r 2 , t) — p(ri, i)/3(r 2 , t)g(ri, r 2 ), 
which can be seen by integrating out the momentum dependence. 

The simplest case is that in which both T := 7I and D := Dol are proportional 
to the identity matrix, when we have the standard definition D = kgT/(m'j). In 
this case it is easy to check that the two formulations agree (see also [3j). This is 
unsurprising as both the difficulties and interest in this analysis lie with the non- 
uniform terms in the friction tensor. 

To demonstrate that the two formulations differ in general, we consider the simple 
example used in Corollary |4. 2 1 In addition, we assume the existence of a parameter A, 
as described above. Then, by the block diagonal form of T, D is also block diagonal 
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with blocks D (l — J^e^i Zi(rj,r^)), i.e. in the notation of [54] wn = — Zi. The result 
to compare with ( |4.2| is (see [54] ) 



9 T p(r,r) = A)V r • [V r p(r,r) + ^p(r, r)V r Fi(r, t) 
+ T^J dry 2 )(r,r',r)V r V/ 2 (r,r') 

-y dr'Zi^rO^^^r'.^ + ^Vr^Cr.^ + ^^rOl^Cr.r'.r) 
+ / drV 3) (r,r',r",r)V r y 2 (r,r") 



Using the approximate two-body form of D in (4.2) gives 



<9 T p(r,r) = D V r • [v r p(r,r) + j^p(r, r)V r ^i(r, r) 
+ drV 2 )(r,r',r)V r T/ 2 (r,r') 

- y dr'ZiCr, r') (p(r', r).g(r, r')V r p(r, r) + ^p( 2 ) (r, r', rJV^^r, r) 
+ ^(r',t). 9 (r,r') / dr'V 2 >(r, r", r)V r V 2 (r, r" 



and it is clear that the two formulations are not, in general, equivalent. See Figure 
|1 . 1| for a diagrammatic representation of the difference in the formalisms. 

Interestingly, despite their obvious differences, in the overdamped limit both for- 
mulations are accurate to 0(e 2 ). There does not seem to be any mathematical or 
physical justification to say that one of them is 'more correct' than the other. How- 
ever, these differences make it clear that the two processes, (i) adiabatically elim- 
inating the fast momentum variable and (ii) integrating over all but one particle's 
coordinates, do not commute. It is worth noting that the need for knowledge of p^ 
in the first case stems from the explicit coupling of the two-body diffusion tensor 
and potential in the iV-body Smoluchowski equation. In contrast, the potential and 
friction tensor are not explicitly coupled in the TV-body Kramers equation, and thus 
only // 2 ) is required. This is a partial explanation of why the resulting equations must 
be different. It would be interesting to perform numerical studies to see if one can 
quantify the differences, i.e. if one can determine the magnitude of the difference in 
the 0(e 2 ) terms. 

We close by stating a result which is most useful when reducing from a phase-space 
dynamical density functional theory to one in only position space: 

Corollary 4.3. Terms proportional to Hermite polynomials of order 2 and 
higher in p enter /(r, p, t) at most with order e 2 . 

5. Conclusions and open problems. Our main result is that, for suitable 
two-body potentials and friction tensors, and using the Enskog approximation in the 



limit of small e, the leading-order solution to (2.5) is given by Theorem 4.1 This is 



a novel Smoluchowski-type equation with a new definition of the one-body diffusion 
tensor. In addition, the Hilbert expansion studied in Section [3] allows us to show 
rigorously that a term typically neglected by heuristic arguments in the derivation of 



DDFT [21 [3] is indeed negligible in the overdamped limit; see Corollary 4.3 However, 
these results have only been shown to hold when the initial condition is independent 
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of e, along with assuming that g is independent of e and p. We now discuss how 
removing these assumptions should be tackled. 

The assumption that the initial condition is independent of e was made for con- 
venience as it allows for analytical progress. In general, however, the question of 
how the initial condition for Kramers equation is related to the correct corresponding 
initial condition for the Smoluchowski equation needs to be addressed. As with the 
Boltzmann equation (cf. [57]) and also noted in the discussion following (3.1), we 
would expect a boundary layer in time (much shorter than the macroscopic timescale 
discussed in this work), over which the given initial condition is attracted to one with 
Gaussian momentum distribution. However, this introduces additional complications, 
that can be studied by modifying appropriately the Hilbert expansion, introducing 
terms that account for the boundary layer and decay exponentially in time, see e.g. 



[5]. We note that even if the non- negativity preserving assumption of Theorem 4.1 



did not hold, then the evolution equations given would be accurate to 0(e). Whilst it 
would be ideal to have a proof that the evolutions in Section [3~4] do preserve the non- 
negativity of ip 1 as mentioned previously, this leads to significant additional technical 
difficulties and may well require further assumptions on the potentials and hydrody- 
namic interactions. 

In addition, as mentioned in Section [T] according to DFT, g is also a functional 
of p. If p is independent of e (i.e. depends only on /o) then the analysis is unaffected. 
This is the case if the initial condition is independent of e and the terms tpi (the part 
of ft in the null space of C) are uniformly zero for all time and all i. As mentioned 
earlier, such a result would rely on the dissipativity of the determining equations, or 
equivalently on the equations being non-negativity preserving. 

If p{v,t) depends on e, then the nonlinear operators are no longer quadratic in / 
as p depends on higher-order, e-dependent parts of /. Furthermore, we do not know 
the precise dependence of g on p. However, if we expand g as a power series in e, 
9 = .9o + «?i + • ■ ■ , then the equations in e~ 2 and e -1 change only by replacing g 
with go- This is because there are no extra terms in the e -2 equation, and only the 
extra term A/"(/o, fo,gi) = enters the e _1 equation. The main point is that g\ enters 
only through non-linear terms, namely by the addition of the terms A/"o(/i, /o, 5i), 
A<o(/o, fx, <7i) and Afi(fo, fa, 9i) to the right hand side of (3.5) (where we have now 
shown the explicit dependence of the non-linear terms on the gj). Thus the conclusion 
that the dynamics of /o are governed by the solvability condition J dp(Cifi-dtfo) = 
still holds. 

The first difference comes when determining / 2 , or more precisely -P1/2, which 
gains additional g\ -dependent terms. The evolution equation for tp looks superficially 
similar, but results in a new definition of a.2 and hence also of &2. 

Allowing g to depend (symmetrically) on pi and P2 introduces many additional 
complications in the analysis. If the p-dependence is introduced at leading order 
it significantly changes the analysis of the non-linear terms. Whilst, by (2.2), it still 
holds that J dr'dp'/W (r, p, r)g(r, r', p, p') = N— 1, we actually require an expr essio n 
for J dp'/W(r, p,r)g(r, r',p, p'). For example, the A/"i(/o,/o) term in Lemma 
significantly more complicated. 

Ideally, one would like to prove a result analogous to Theorem 3.1 of [27], which 
states that, under suitable assumptions, in a suitable norm, and for a fixed macro- 
scopic time period (0, to], the solution to the Boltzmann equation is 0(e) close to the 
local Maxwellian whose parameters vary according to the hydrodynamic equations. 
Note that in the case where to —¥ 00, the constant in the 0(e) bound may diverge. 
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3.9 



is 



In order to prove such result, one must truncate the Hilbert series at a finite order 
and add a remainder term. One then determines bounds on each of these terms, which 
require sufficiently good estimates on the collision term (the hydrodynamic interac- 
tions and V2 in our case). This truncation is necessary as the Hilbert expansion does 
not converge uniformly in the small parameter. Since such estimates on the collision 
operator depend on its precise form (in particular, it is assumed that the kernel of 
the Boltzmann collision operator has finite range; not true for hydrodynamic interac- 
tions), and a specific choice of norms, we have restricted our analysis to determining 
the leading order terms in such an expansion. 

We close by discussing some open problems. The first area concerns confined 
fluids and the effects of boundaries. Although the external potential V\ may be used 
to model boundaries which are impermeable to the colloid particles but permeable 
to the fluid, a truly confined fluid cannot be modelled in this way. Extension to 
such systems would require a treatment of the hydrodynamic interactions caused by 
the boundaries. Such effects break the symmetry of the bath, as well as changing 
the mobility of the colloid particles near the boundaries. Additional complications 
would result from the presence of heterogeneities at boundaries, which is indeed the 
case in practice. Heterogeneous boundaries, either chemical or topographical, can 
have a significant effect on the behavior of fluids both at both the micro-scale (e.g. 
they can influence the thickness of the wetting layer in the immediate vicinity of the 
boundary and corresponding wetting transitions) and macro-scale (they can affect the 
shape of the gas-liquid interface away from the boundaries) [SSJ H2J 03 ISH1 (Ml US] ■ It 
would also be of interest to study mixtures of colloid particles, e.g. a system with 
two types of particle which differ in their sizes, masses, or interparticle potentials 
V-2. As mentioned above, a full treatment of the problem would involve analysis 
of boundary layer effects, including how the initial condition for the Smoluchowski 
equation should be determined by that for the Kramers equation. These and related 
issues are currently under investigation. 
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Appendix A. Proofs of some lemmas of Section [3} 



A.l. Proof of Lemma 3.1, Proof. Let w = (v(ri, pi, i), . . . v(rjv,p^, t)J 
Since r(r w ) is positive-definite, we have, for some 8 > 0, w • Tw > 5w ■ w. Hence 



N 

W • Tw = ^ v ( r iiPi>*) ' r ij( rJV ) v ( r i>P^) 

N 

= X^ V ( ri ' Pl '^ ' 1 + ^2 Z l( V i> r j) v ( r »)Pi)*) 

i=l jjti 

N 

+ ^Z V ( ri ' Pi '^ ' Z 2( r »> r i) v ( r i>Pi>*) ^ S Y1 l V ( r i'Pi'*)P 

i^tj i=l 
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T 



Since f^ N>> (r N ,p N ,t) is non-negative, and positive on a set of non-zero measure (as 
by definition /W > and / dr N dp N f W = jv), we have 



N 



=1 



•53v(ri,pi,t) • Z 2 (r i ,r i )v(r 3 -,p J -,t) 



AT 

2 — 1 

By the symmetry of f^ N \ interchanging dummy variables of integration gives 
J dr%"/W(r^p",i)(jVv(r 1)Pl ,i). [l + (iV - f )Z 1 (r 1 , r 2 )] v(n, Pl , 

+ N(N - l)v(ri, pi, t) • Z 2 (r l5 r 2 )v(r 2 , p 2 , t) 
><57V| dr N dp K /( N )(r N ,p N ,t)|v(r llPll ()| 2 
Using (2.2 1 for the cases with n = 1 and n = 2, i.e. 

/( 2 )(r 1)Pl ,r 2 ,p 2) i)=iV(iV-l) | dr JV - 2 dp A, - 2 / W (r A ', P JV ,t)and 
/«(r 1)Pl ,i)=JV | dr^dp^/^V.pV) 

gives 

y dridpidr 2 dp 2 / (2) (ri,pi,r 2 ,p 2 ,t)^v(ri,pi,t) • + Zi(r 1; r 2 ) v(ri,pi,i) 

+ v(ri,pi,t) • Z 2 (ri,r 2 )v(r 2 ,p 2 ,t)^ 

> 6 J dridpi/(ri,pi,i)|v(ri,pi,f)| 2 . 

Inserting the definition /( 2 )(n,pi,r 2 ,p 2 ,t) = g(r lf r 2 )f(r u pi, i)/(r 2 , p 2 , *) and re- 
naming the dummy variables gives the result. 

The fact that / may be chosen as is tr ivial. To see that the result holds when / 
is replaced by /o, we insert the expansion (3.1 ) and note / dr'dp'g(r, r')/^ 1 '^', p' , t) = 
N — 1 holds for all e, in particular for e = 0. □ 



A. 2. Proof of Lemma |3.6[ Proof. We consider each of the three op erat ors 
in C individually, starting with Cq. For arbitrary /, /, and using Corollary 
particular that V p (/ _1 ) = p/q 1 



3.3 



</,A)/) /o -i = / drdp/ - 1 (r,p,<)/(r, P ,t)V p -[(p + V p )/(r,p,t)] 

= - J drdp/,7 1 (r, p, t)(p + V p )/(r, p, t) ■ [(p + V p )/(r, p, t)] 
drdp/o" 1 ^, p, i)/(r, p, *)V P • [(p + V p )/(r, p, t)], 
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where the second and third lines both follow via integration by parts. Hence Cq is 
self- adjoint. 



For Af (f , /), by Lemma 



3.4 



we have 



</, Wo, /))/-* = / drdp/ - 1 (r,p,t)/(r,p s t) 



x V p - [ / dr'dp' 3 (r,r')/ (r',p',t)Z 1 (r,r') x (p + V p )/(r, p, t) 

= - J drdp/Q-^r, p, i)(p + V p )/(r, p, t) 

■ J dr'dp'. 9 (r,r')/o(r',p',t)Z 1 (r,r') x (p + V p )/(r, p, t) 

= - J drdpdr>7 -\r,p,t). 9 (ry)/ (r',p',i) 

x (p + V p )/(r, p, t) ■ Zx(r, r')(p + V p )/(r, p, t), 

where we have used integration by parts and Fubini's theorem. Now, since Zi is a 
symmetric matrix, / and / can be interchanged and the argument reversed, showing 
that A/q(/o, /) is self-adjoint. 



It remains to calculate the adjoint of Aq(/, Jo)- Using Lemma 



3.4 



gives 



(/,W,/o))f„- 1= - / drdp/ -V,p,i)/(r,p,f)/o(r,p,t) 

dr'dp' 3 (r, r')Z 2 (r, r')(p' + V P 0/>', p', t) ■ p 
drdpdr'dp' ff (r, r')/(r, p, t)f(r', p', t)p ■ Z 2 (r, r')p', 



where we have used the divergence theorem, Fubini's theorem and the identity (for 
symmetric matrices) Z 2 p' ■ p = p ■ Z 2 p'. Since this final term, along with the rest of 
the integral is symmetric under interchanging the pairs of dummy variables (r, p) ■<-> 
(r',p') we see that Ao(/,/o) is also self-adjoint. The overall result now follows from 
linearity of the integral, and hence of the adjoint. □ 



A. 3. Proof of Lemma 3.7, Proo f. W e prove the equivalent statement (for 

self adjoint operators, as C is by Lemma 3.6) that there exists an orthonormal basis 

of L 2 i such that ££» = A,£j, Xj € K such that lim,_ J . 00 |A, | = oo iTS Theorem 

J Jo 

11.3.13]. We make use of Lemma 3.1 which allows us to compare the eigenvalues of 

L to those where Zj = 0, and the fact that the eigenf unctions and eigenvalues of the 

resulting operator can be constructed explicitly. We note that the P„, a form a basis 

of L 2 (R 3 , e-IPl 2 / 2 ) and so the f unctions e -l p l 2 / 2 P„, a form a basis of L 2 (R 3 , e lF»l 2 /2 ). 

First note, by Lemma |3.12| that the spaces 

Span{e" |p|2/2 P„, a (p)|n fixed, |a| = n} 

(where the coefficients may be functions of r, t) are invariant under C. Thus all 
eigenfunctions may be written in the form 



Vy, (r, p, t) = £a,j(r, t) e^ 2 ' 2 P„, a (p), 

a 
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(A.l) 



where j = 1, . . . , T(n + 1), with T(n) the n-th triangular number (which corresponds 
the the numb er of solutions to a% + a 2 + a 3 — n — 1). Also, as noted in the proof of 



Lemma 3.12 7V(/, /o) contributes only for n = 1, and as such we may ignore it when 
calculating the eigenvalues. It therefore suffices to consider the eigenvalues of 

Cf(r, p, t) = V p • [ J dr'dp'.g(r, r')/o(r', p', t)(j^l + Z x (r, r')) (p + V p )/(r, p, t) 

=: V p -Z(r,f)(p + V p )/(r,p,t). 
Suppose —Cipn.j = \ n ,jipn,ji then 

X nJ J drdp/^ 1 (r, p, t)\ip nJ (r, p, i) | 2 

= - J drdp/ ( 7 1 (r,p,i)V'„,j(r,p,t)Vp • Z(r, i)(p + Vp)^n,j(r, p, i) 
= / drdp/- 1 (r, p, t) (p + V p )Vnj (r, P, *) ' Z(r, t) (p + V p )V»j (r, P, *) 

= y drdpdr / dp7o- 1 (r ! P,t)5(r,r / )/o(r / ! p',t) 

x (p + V p )^(r,p,t) • ( ] ^ T l + Z 1 (r,r / ))(p + Vp)V„j(r,p,t) 
= /" drdpdr'dp'/o^p^)^^,?',^^^) 

x v(r, p, t) ■ (j^l + Zx(r, r')) v(r, p, t), 

where the second equality follows via integration by parts and that / ( ^~ 1 (r,p,<) = 
Pq 1 (r,i)Zexp(|p| 2 /2), and we denote v(r,p,i) = (p + V p )V'n,j(i',P,t)/o" 1 ( r »P>*)- 
Now note that Lemma [3 . 1 1 holds when Z 2 is set to zero since it requires only that T 
is positive definite with the correct symmetry. Since Tn is a principal minor of T, 
it is positive definite, and by symmetry so are all Tjj. It therefore follows that the 
block diagonal matrix with entries Tjj is also positive definite, with the same required 
symmetry as T and we have 

y drdp/,f 1 (r,p,t)|V>,,,j(r,p,t)| 2 

><5 y drdp/ (r,p,t)|v(r,p,t)| 2 ^ J drdp/^r, p, t)|(p + V p )VVj <(r, P, t)\ 2 

= -$[ drdp/ ( 7 1 (r,p,i)V', i j(r,p,t)Vp • (p + V p )-0„j(r,p,t). 
We now compute V p • (p + ^7 p )tp n j(r,p,t): 

3 

Vp ■ (p + V p ) e-!Pl 2 /2 Pn a(p) = d Pj ( Pj + d Pj ) e-IPl 2 / 2 P„, a (p) 

3=1 

= E ^ ( e ' P ' 2/2 d K F ».^P)) = E ^ ( e ' P ' 2/2 ",/',a , (P)) 

3=1 3=1 

= £ elPl 2 / 2 a,-(- Fj +a Pj )P„, a _e 3 .(p) = -elPl 2 / 2 P„ a (p)^a, 

3=1 3=1 

= -nclPl 2 / 2 P„, a (p). 
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The required identities for operators on P„ !a follow from its product form and the 
equivalent 1-dimensional identities. Note Bj is the jth unit vector. Thus we have 



J drdp/ 1 (r,p > t)|^nj(r,p,t)| a >nS J drdp/ l (r, p, t)\i/> nJ (r, p, t)f 



and the result follows. □ 



A. 4. Proof of Lemma 3.12, Proof. From (2.6a) and Lemma 3.4 for CqJ and 



■M)(/o> /), it is sufficient to consider a general operator 

£ := V p ■ Z(r)(p + V p ) = £ ^0 Pi fe + Pj ). 

We have, using the standard identity d x H n (x) = nH n _i(x), 

(p + d p )[H a (p)e-P 2 ^} = e~v 2 ' 2 d p H a {p) = e^ 2 / 2 aH a ^(p) (A.2) 

and using H n+ i(x) = xH n (x) — d x H n (x), we find 

d p [H a (p)e'P 2 / 2 } = a-? 2 / 2 (-p + d p )H a (p) = -H a+1 (p)e-P 2 / 2 

It is therefore clear that C preserves |a| = n, with the possibility of the new coefficients 
all being zero. 



It remains to consider Ao(/, /o), which by Lemma 3.4 is given by 



W, Jo) = -^4r/o(r, p, t) J dr'dp' 5 (r, r')Z 2 (r, r')(p' + V P 0/(r', p', t) ■ p 



Note that P = 1. Using |A.2| , f a = Z" 1 e~ p /( 2mk BT) po ( T> t), and that the P n ^ are 
orthogonal, it is clear that the integral gives zero for any terms not proportional to 
Pi = Hn, and in this case returns something of the form a(r) ■ p. 

Hence C preserves n and it remains to show that £ X}|a|=n a n,a{ r ,p)P n ,a(p) = if 
and only if a n ^ a = for all a. This follows from the null space of C being e - ' p ' 2 / 2 <p(r, t) 



see Lemma 3.8 1, and thus containing only P , and the orthogonality of the P n ,&- □ 
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